      SUBROUTINE ACRSSB(A,B,X,N)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     SUBROUTINE ACRSSB                                                C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DIMENSION A(N),B(N),X(N)
C
      DO 1 I=1,N
      J=I+1
      IF (I.EQ.N) J=1
      K=J+1
      IF (J.EQ.N) K=1
    1 X(I)=A(J)*B(K)-A(K)*B(J)
C
      RETURN
      END
      FUNCTION ADOTB(A,B,N)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     FUNCTION ADOTB                                                   C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DIMENSION A(N),B(N)
      S=0
      DO 1 I=1,N
    1 S=S+A(I)*B(I)
      ADOTB=S
      END
      FUNCTION ALENGT(P,N)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     FUNCTION ALENGT                                                  C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DIMENSION P(N)
      D=0.E0
      DO 1 I=1,N
    1 D=D+P(I)**2
      ALENGT=SQRT(D)
      RETURN
      END
      FUNCTION ANGLE(C1,C2,C3,X,Y,N)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     FUNCTION ANGLE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DIMENSION C1(N),C2(N),C3(N),X(N),Y(N)
C
      DO 10 I=1,N
      X(I)=C1(I)-C2(I)
   10 Y(I)=C3(I)-C2(I)
      ANGLE=ACOS(ADOTB(X,Y,N)/ALENGT(X,N)/ALENGT(Y,N))
C
      RETURN
      END
      FUNCTION DIHEDR(C1,C2,C3,C4,AX,X,Y,Z,N)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     FUNCTION DIHEDR                                                 C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      AXL=DSTNC2(C2,C3,AX,N)
C
      CALL DIF(C1,C2,X,N)
      A=ADOTB(X,AX,N)/AXL
      CALL SCALE(AX,A,Z,N)
      CALL DIF(X,Z,X,N)
C
      CALL DIF(C4,C3,Y,N)
      A=ADOTB(Y,AX,N)/AXL
      CALL SCALE(AX,A,Z,N)
      CALL DIF(Y,Z,Y,N)
C
      CALL ACRSSB(X,Y,Z,N)
      D1=ADOTB(X,Y,N)/ALENGT(X,N)/ALENGT(Y,N)
      IF (ABS(D1) .GT. 1.0) D1=NINT(D1)
      A=ACOS(D1)
      DIHEDR=-SIGN(A,ADOTB(AX,Z,N))
C
      RETURN
      END
      FUNCTION DIST(P,Q)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     FUNCTION DIST                                                    C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DIMENSION P(3),Q(3)
      D=0
      DO 1 I=1,3
    1 D=D+(P(I)-Q(I))**2
      DIST=SQRT(D)
      RETURN
      END
      FUNCTION DIST2(P,Q,N)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     FUNCTION DIST2                                                   C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DIMENSION P(N),Q(N)
      D=0
      DO 1 I=1,N
    1 D=D+(P(I)-Q(I))**2
      DIST2=D
      RETURN
      END
      FUNCTION DSTNCE(P,Q,X)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     FUNCTION DSTNCE                                                  C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DIMENSION P(3),Q(3),X(3)
      D=0.E0
      DO 1 I=1,3
      X(I)=P(I)-Q(I)
    1 D=D+X(I)**2
      DSTNCE=SQRT(D)
      END
      FUNCTION DSTNC2(P,Q,X,N)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     FUNCTION DSTNC2                                                  C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DIMENSION P(N),Q(N),X(N)
      D=0.E0
      DO 1 I=1,N
      X(I)=P(I)-Q(I)
    1 D=D+X(I)**2
      DSTNC2=D
      END
      SUBROUTINE DIF(P,Q,X,N)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     SUBROUTINE DIF                                                   C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DIMENSION P(N),Q(N),X(N)
C
      DO 1 I=1,N
    1 X(I)=P(I)-Q(I)
C
      RETURN
      END
      SUBROUTINE SUM(P,Q,X,N)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     SUBROUTINE SUM                                                   C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DIMENSION P(N),Q(N),X(N)
C
      DO 1 I=1,N
    1 X(I)=P(I)+Q(I)
C
      RETURN
      END
      SUBROUTINE INDEXX(N,ARRIN,INDX)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     Indexes the array ARRIN of length N, i.e. outputs the array INDX
C     such that ARRIN(INDX(J)) is in ascending order for J=1,2,...,N.
C     The input quantities N and ARRIN are not changed.
C
C     Source: Numerical Recipes 
C     by W.H.Press, B.P.Flannery, S.A.Teukolsky and W.T.Vetterling
C     Cambridge University Press, 1986.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IMPLICIT INTEGER (A-Z)
      DIMENSION ARRIN(N),INDX(N)
      DO 11 J=1,N
        INDX(J)=J
11    CONTINUE
C----- THE FOLLOWING ADDED TO ENSURE VALIDITY FOR N=1
      IF (N.EQ.1) RETURN
      L=N/2+1
      IR=N
10    CONTINUE
        IF(L.GT.1)THEN
          L=L-1
          INDXT=INDX(L)
          Q=ARRIN(INDXT)
        ELSE
          INDXT=INDX(IR)
          Q=ARRIN(INDXT)
          INDX(IR)=INDX(1)
          IR=IR-1
          IF(IR.EQ.1)THEN
            INDX(1)=INDXT
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF(J.LE.IR)THEN
          IF(J.LT.IR)THEN
            IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
          ENDIF
          IF(Q.LT.ARRIN(INDX(J)))THEN
            INDX(I)=INDX(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GO TO 20
        ENDIF
        INDX(I)=INDXT
      GO TO 10
      END
      SUBROUTINE INDEXR(N,ARRIN,INDX)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     Version of INDEXX with REAL array ARRIN.
C
C     Indexes the array ARRIN of length N, i.e. outputs the array INDX
C     such that ARRIN(INDX(J)) is in ascending order for J=1,2,...,N.
C     The input quantities N and ARRIN are not changed.
C
C     Source: Numerical Recipes 
C     by W.H.Press, B.P.Flannery, S.A.Teukolsky and W.T.Vetterling
C     Cambridge University Press, 1986.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IMPLICIT INTEGER (A-Z)
      REAL ARRIN(N),Q
      DIMENSION INDX(N)
      DO 11 J=1,N
        INDX(J)=J
11    CONTINUE
C----- THE FOLLOWING ADDED TO ENSURE VALIDITY FOR N=1
      IF (N.EQ.1) RETURN
      L=N/2+1
      IR=N
10    CONTINUE
        IF(L.GT.1)THEN
          L=L-1
          INDXT=INDX(L)
          Q=ARRIN(INDXT)
        ELSE
          INDXT=INDX(IR)
          Q=ARRIN(INDXT)
          INDX(IR)=INDX(1)
          IR=IR-1
          IF(IR.EQ.1)THEN
            INDX(1)=INDXT
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF(J.LE.IR)THEN
          IF(J.LT.IR)THEN
            IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
          ENDIF
          IF(Q.LT.ARRIN(INDX(J)))THEN
            INDX(I)=INDX(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GO TO 20
        ENDIF
        INDX(I)=INDXT
      GO TO 10
      END
      SUBROUTINE MATSUM(A,B,AFACT,BFACT,C,N)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     SUBROUTINE MATSUM                                                C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DIMENSION  A(N),B(N),C(N)
      DO 1 I=1,N
    1 C(I)=AFACT*A(I)+BFACT*B(I)
      RETURN
      END
      SUBROUTINE MINV (A,N,D,L,M)
C
CCCCCC W.F. VAN GUNSTEREN, CAMBRIDGE, APR. 1979 CCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     SUBROUTINE MINV (A,N,D,L,M)                                      C
C                                                                      C
COMMENT   MINV INVERTS THE N*N MATRIX A, USING THE STANDARD GAUSS-     C
C     JORDAN METHOD. THE DETERMINANT IS ALSO CALCULATED. IT HAS BEEN   C
C     COPIED FROM THE IBM SCIENTIFIC SUBROUTINE PACKAGE.               C
C                                                                      C
C     A(1..N,1..N) = MATRIX TO BE INVERTED; DELIVERED WITH THE INVERSE C
C     N = ORDER OF MATRIX A                                            C
C     D = DELIVERED WITH THE DETERMINANT OF A                          C
C     L,M(1..N) = DUMMY ARRAYS                                         C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DIMENSION A(1),L(1),M(1)
C
C*****SEARCH FOR LARGEST ELEMENT
      D=1.E0
      NK=-N
      DO 80 K=1,N
      NK=NK+N
      L(K)=K
      M(K)=K
      KK=NK+K
      BIGA=A(KK)
      DO 20 J=K,N
      IZ=N*(J-1)
      DO 20 I=K,N
      IJ=IZ+I
   10 IF ( ABS(BIGA)- ABS(A(IJ))) 15,20,20
   15 BIGA=A(IJ)
      L(K)=I
      M(K)=J
   20 CONTINUE
C
C*****INTERCHANGE ROWS
      J=L(K)
      IF (J-K) 35,35,25
   25 KI=K-N
      DO 30 I=1,N
      KI=KI+N
      HOLD=-A(KI)
      JI=KI-K+J
      A(KI)=A(JI)
   30 A(JI)=HOLD
C
C*****INTERCHANGE COLUMNS
   35 I=M(K)
      IF (I-K) 45,45,38
   38 JP=N*(I-1)
      DO 40 J=1,N
      JK=NK+J
      JI=JP+J
      HOLD=-A(JK)
      A(JK)=A(JI)
   40 A(JI)=HOLD
C
C*****DIVIDE COLUMN BY MINUS PIVOT
   45 IF (BIGA) 48,46,48
   46 D=0.E0
      RETURN
   48 DO 55 I=1,N
      IF (I-K) 50,55,50
   50 IK=NK+I
      A(IK)=A(IK)/(-BIGA)
   55 CONTINUE
C
C*****REDUCE MATRIX
      DO 65 I=1,N
      IK=NK+I
      HOLD=A(IK)
      IJ=I-N
      DO 65 J=1,N
      IJ=IJ+N
      IF (I-K) 60,65,60
   60 IF (J-K) 62,65,62
   62 KJ=IJ-I+K
      A(IJ)=HOLD*A(KJ)+A(IJ)
   65 CONTINUE
C
C*****DIVIDE ROW BY PIVOT
      KJ=K-N
      DO 75 J=1,N
      KJ=KJ+N
      IF (J-K) 70,75,70
   70 A(KJ)=A(KJ)/BIGA
   75 CONTINUE
C
C*****PRODUCT OF PIVOTS
      D=D*BIGA
C
C*****REPLACE PIVOT BY RECIPROCAL
      A(KK)=1.E0/BIGA
   80 CONTINUE
C
C*****FINAL ROW AND COLUMN INTERCHANGE
      K=N
  100 K=(K-1)
      IF (K) 150,150,105
  105 I=L(K)
      IF (I-K) 120,120,108
  108 JQ=N*(K-1)
      JR=N*(I-1)
      DO 110 J=1,N
      JK=JQ+J
      HOLD=A(JK)
      JI=JR+J
      A(JK)=-A(JI)
  110 A(JI)=HOLD
  120 J=M(K)
      IF (J-K) 100,100,125
  125 KI=K-N
      DO 130 I=1,N
      KI=KI+N
      HOLD=A(KI)
      JI=KI-K+J
      A(KI)=-A(JI)
  130 A(JI)=HOLD
      GOTO 100
  150 RETURN
      END
      SUBROUTINE NORM(P,N)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     SUBROUTINE NORM                                                  C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DIMENSION P(N)
C
      D=0.E0
      DO 10 I=1,N
   10 D=D+P(I)**2
      D=1.E0/SQRT(D)
      DO 20 I=1,N
   20 P(I)=P(I)*D
C
      RETURN
      END
      SUBROUTINE RANDOM (RAND,IG)
C
CCCCCC R. GEURTSEN, GRONINGEN, WFVG, JULY 1987 CCCCCCCCCCCCCCCCCCCCCCCCC
C      MODIFIED BY: TON RULLMANN, UTRECHT, FEBRUARY 1990               C
C                                                                      C
C     SUBROUTINE RANDOM (RAND,IG)                                      C
C                                                                      C
COMMENT   RANDOM GENERATES A RANDOM NUMBER RAND, USING A LINEAR        C
C     CONGRUENTIAL METHOD. THE RECURSION FORMULA                       C
C                                                                      C
C         IRAND = MOD(IRAND * B + 1, A)                                C
C                                                                      C
C     IS USED WITH  B = 31415821  AND  A = 100000000. THE LAST DIGIT   C
C     FROM THE RANDOM INTEGER IRAND IS CHOPPED OF, AND THE NUMBER      C
C     IS SCALED TO A REAL VALUE RAND BETWEEN 0 AND 1, INCLUDING 0 BUT  C
C     EXCLUDING 1.                                                     C
C                                                                      C
C     RAND = DELIVERED WITH RANDOM NUMBER BETWEEN 0 AND 1              C
C     IG = RANDOM NUMBER GENERATOR SEED, IS DELIVERED WITH RANDOM      C
C          INTEGER                                                     C
C                                                                      C
C MODIFIED: ON SUBSEQUENT ENTRIES IRAND IS AGAIN CALCULATED FROM IG,   C
C     JUST AS ON THE FIRST ENTRY. FOR A CONTINUED SERIES THE USER MUST C
C     PROVIDE THE OUTPUT VALUE OF IRAND AS INPUT TO THE NEXT CALL.     C
C MODIFIED: PRINT ERROR MESSAGE IF OUT OF RANGE (INCLUDE FILE NEEDED)  C
C                                                                      C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      INCLUDE 'STANDFIL.INC'
      PARAMETER (M=100000000, M1=10000, MULT=31415821)
      DATA       IRAND /0/, NEW /0/
C
C
C*****TR: NEXT LINE COMMENTED, SO IRAND IS RECALCULATED FROM IG
C      IF (NEW .NE. 0) GOTO 7
      NEW = 1
      IRAND = MOD (IABS(IG),M)
    7 CONTINUE
C
C*****MULTIPLY IRAND BY MULT, BUT TAKE INTO ACCOUNT THAT OVERFLOW MUST
C*****BE DISCARDED, AND DO NOT GENERATE AN ERROR.
C
      IRANDH = IRAND / M1
      IRANDL = MOD(IRAND, M1)
      MULTH = MULT / M1
      MULTL = MOD(MULT, M1)
C
      IRAND = MOD(IRANDH*MULTL + IRANDL*MULTH, M1) * M1 + IRANDL*MULTL
      IRAND = MOD(IRAND + 1, M)
C
C*****CONVERT IRAND TO A REAL RANDOM NUMBER BETWEEN 0 AND 1.
C
      R = REAL(IRAND / 10) * 10 / REAL(M)
      IF ((R .LT. 0.E0) .OR. (R .GE. 1.E0)) THEN
         R = 0.E0
         WRITE (ITXERR,*) 
     +     ' *** RANDOM: RANDOM NUMBER OUT OF RANGE, SET TO 0.'
      ENDIF
      RAND = R
      IG = IRAND
C
      RETURN
      END
      FUNCTION RMS(X2,X)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     FUNCTION RMS                                                     C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      Y=X2-X*X
      IF (Y.LT.0.E0) Y=0.E0
      RMS=SQRT(Y)
      RETURN
      END
      SUBROUTINE ROTATE(A,N,INDEX,X,Y,Z)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     SUBROUTINE ROTATE                                                C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      CHARACTER*(*) INDEX
      DIMENSION A(3,N),P(3,3),Q(3,3),R(3,3)
C
   10 FORMAT(1X,A10,' ROTATION ANGLE FOR X,Y,Z AXIS (REPLACED BY',
     + ' RANDOM NUMBER IF -1.)')
      IF (INDEX.NE.' '.AND.X.EQ.0..AND.Y.EQ.0..AND.Z.EQ.0.) THEN
        PRINT 10,INDEX
        READ *,X,Y,Z
        CALL RANGEN
        IF(X.EQ.-1.) X=RANGET()*3.14
        IF(Y.EQ.-1.) Y=RANGET()*3.14
        IF(Z.EQ.-1.) Z=RANGET()*3.14
      END IF
      CALL ZERO(P,9)
      CALL ZERO(Q,9)
      P(1,1)=1.
      P(3,3)=COS(X)
      P(2,2)=P(3,3)
      P(2,3)=SIN(X)
      P(3,2)=-P(2,3)
      Q(2,2)=1.
      Q(3,3)=COS(Y)
      Q(1,1)=Q(3,3)
      Q(1,3)=SIN(Y)
      Q(3,1)=-Q(1,3)
      DO 100 I=1,3
      DO 100 J=1,3
      R(I,J)=0.
      DO 90 K=1,3
   90 R(I,J)=R(I,J)+P(I,K)*Q(K,J)
  100 CONTINUE
      CALL ZERO(P,9)
      P(3,3)=1.
      P(2,2)=COS(Z)
      P(1,1)=P(2,2)
      P(1,2)=SIN(Z)
      P(2,1)=-P(1,2)
      DO 150 I=1,3
      DO 150 J=1,3
      Q(I,J)=0.
      DO 140 K=1,3
  140 Q(I,J)=Q(I,J)+P(I,K)*R(K,J)
  150 CONTINUE
C     PRINT 160,Q
  160 FORMAT(' ROTATION MATRIX:'/3F10.4)
      DO 200 L=1,N
      DO 180 I=1,3
      P(I,1)=0.
      DO 170 J=1,3
  170 P(I,1)=P(I,1)+Q(I,J)*A(J,L)
  180 CONTINUE
      DO 190 I=1,3
  190 A(I,L)=P(I,1)
  200 CONTINUE
      RETURN
      END
      SUBROUTINE SCALE(P,SC,X,N)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     SUBROUTINE SCALE                                                 C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DIMENSION P(N),X(N)
C
      DO 1 I=1,N
    1 X(I)=P(I)*SC
C
      RETURN
      END
      FUNCTION SPUR(A,B,N)
C********************************************************************
C     FUNCTION SPUR
C********************************************************************
C
C     SPUR COMPUTES TRACE OF MATRIX PRODUCT A*B
C     OF SYMMETRIC MATRICES, THE LOWER TRIANGLES OF WHICH ARE GIVEN
C     AS VECTORS IN A AND B, RESPECTIVELY.
C
      DIMENSION A(*),B(*)
      K=0
      SPUR=0.
      DO 10 I=1,N
      DO 10 J=1,I
      K=K+1
      S=A(K)*B(K)
      SPUR=SPUR+S+S
      IF (I.EQ.J) SPUR=SPUR-S
   10 CONTINUE
      RETURN
      END
      SUBROUTINE TRED2(NM, N, A, D, E, Z)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     SUBROUTINE TRED2                                                 C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DIMENSION A(NM,N), D(N), E(N), Z(NM,N)
      DO 20 I=1,N
C
         DO 10 J=1,I
            Z(I,J) = A(I,J)
10       CONTINUE
20    CONTINUE
C
      IF (N.EQ.1) GO TO 160
C     ********** FOR I=N STEP -1 UNTIL 2 DO -- **********
      DO 150 II=2,N
         I = N + 2 - II
         L = I - 1
         H = 0.0
         SCALE = 0.0
         IF (L.LT.2) GO TO 40
C     ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) **********
         DO 30 K=1,L
            SCALE = SCALE + ABS(Z(I,K))
30       CONTINUE
C
         IF (SCALE.NE.0.0) GO TO 50
40       E(I) = Z(I,L)
         GO TO 140
C
50       DO 60 K=1,L
            Z(I,K) = Z(I,K)/SCALE
            H = H + Z(I,K)*Z(I,K)
60       CONTINUE
C
         F = Z(I,L)
         G = -SIGN(SQRT(H),F)
         E(I) = SCALE*G
         H = H - F*G
         Z(I,L) = F - G
         F = 0.0
C
         DO 100 J=1,L
            Z(J,I) = Z(I,J)/(SCALE*H)
            G = 0.0
C     ********** FORM ELEMENT OF A*U **********
            DO 70 K=1,J
               G = G + Z(J,K)*Z(I,K)
70          CONTINUE
C
            JP1 = J + 1
            IF (L.LT.JP1) GO TO 90
C
            DO 80 K=JP1,L
               G = G + Z(K,J)*Z(I,K)
80          CONTINUE
C     ********** FORM ELEMENT OF P **********
90          E(J) = G/H
            F = F + E(J)*Z(I,J)
100      CONTINUE
C
         HH = F/(H+H)
C     ********** FORM REDUCED A **********
         DO 120 J=1,L
            F = Z(I,J)
            G = E(J) - HH*F
            E(J) = G
C
            DO 110 K=1,J
               Z(J,K) = Z(J,K) - F*E(K) - G*Z(I,K)
110         CONTINUE
120      CONTINUE
C
         DO 130 K=1,L
            Z(I,K) = SCALE*Z(I,K)
130      CONTINUE
C
140      D(I) = H
150   CONTINUE
C
160   D(1) = 0.0
      E(1) = 0.0
C     ********** ACCUMULATION OF TRANSFORMATION MATRICES **********
      DO 220 I=1,N
         L = I - 1
         IF (D(I).EQ.0.0) GO TO 200
C
         DO 190 J=1,L
            G = 0.0
C
            DO 170 K=1,L
               G = G + Z(I,K)*Z(K,J)
170         CONTINUE
C
            DO 180 K=1,L
               Z(K,J) = Z(K,J) - G*Z(K,I)
180         CONTINUE
190      CONTINUE
C
200      D(I) = Z(I,I)
         Z(I,I) = 1.0
         IF (L.LT.1) GO TO 220
C
         DO 210 J=1,L
            Z(I,J) = 0.0
            Z(J,I) = 0.0
210      CONTINUE
C
220   CONTINUE
C
      RETURN
      END
      SUBROUTINE TQL2(NM, N, D, E, Z, IERR)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     SUBROUTINE TQL2                                                  C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DIMENSION D(N), E(N), Z(NM,N)
      REAL MACHEP
C                **********
      MACHEP = 2.**(-22)
C     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
C     ABOVE VALUE SUITABLE FOR VAX-VMS REAL*4 ARITHMETIC.
C
      IERR = 0
      IF (N.EQ.1) GO TO 160
C
      DO 10 I=2,N
         E(I-1) = E(I)
10    CONTINUE
C
      F = 0.0
      B = 0.0
      E(N) = 0.0
C
      DO 110 L=1,N
         J = 0
         H = MACHEP*(ABS(D(L))+ABS(E(L)))
         IF (B.LT.H) B = H
C     ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
         DO 20 M=L,N
            IF (ABS(E(M)).LE.B) GO TO 30
C     ********** E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
C                THROUGH THE BOTTOM OF THE LOOP **********
20       CONTINUE
30       IF (M.EQ.L) GO TO 100
40       IF (J.EQ.30) GO TO 150
         J = J + 1
C     ********** FORM SHIFT **********
         P = (D(L+1)-D(L))/(2.0*E(L))
         R = SQRT(P*P+1.0)
         H = D(L) - E(L)/(P+SIGN(R,P))
C
         DO 50 I=L,N
            D(I) = D(I) - H
50       CONTINUE
C
         F = F + H
C     ********** QL TRANSFORMATION **********
         P = D(M)
         C = 1.0
         S = 0.0
         MML = M - L
C     ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
         DO 90 II=1,MML
            I = M - II
            G = C*E(I)
            H = C*P
            IF (ABS(P).LT.ABS(E(I))) GO TO 60
            C = E(I)/P
            R = SQRT(C*C+1.0)
            E(I+1) = S*P*R
            S = C/R
            C = 1.0/R
            GO TO 70
60          C = P/E(I)
            R = SQRT(C*C+1.0)
            E(I+1) = S*E(I)*R
            S = 1.0/R
            C = C*S
70          P = C*D(I) - S*G
            D(I+1) = H + S*(C*G+S*D(I))
C     ********** FORM VECTOR **********
            DO 80 K=1,N
               H = Z(K,I+1)
               Z(K,I+1) = S*Z(K,I) + C*H
               Z(K,I) = C*Z(K,I) - S*H
80          CONTINUE
C
90       CONTINUE
C
         E(L) = S*P
         D(L) = C*P
         IF (ABS(E(L)).GT.B) GO TO 40
100      D(L) = D(L) + F
110   CONTINUE
C     ********** ORDER EIGENVALUES AND EIGENVECTORS **********
      DO 140 II=2,N
         I = II - 1
         K = I
         P = D(I)
C
         DO 120 J=II,N
            IF (D(J).GE.P) GO TO 120
            K = J
            P = D(J)
120      CONTINUE
C
         IF (K.EQ.I) GO TO 140
         D(K) = D(I)
         D(I) = P
C
         DO 130 J=1,N
            P = Z(J,I)
            Z(J,I) = Z(J,K)
            Z(J,K) = P
130      CONTINUE
C
140   CONTINUE
C
      GO TO 160
C     ********** SET ERROR -- NO CONVERGENCE TO AN
C                EIGENVALUE AFTER 30 ITERATIONS **********
150   IERR = L
160   RETURN
C     ********** LAST CARD OF TQL2 **********
      END
      SUBROUTINE TRPMXV(AMX,VEC,RES,NR,NC)
C********************************************************************
C
C     SUBROUTINE TRPMXV PERFORMS A MATRIX*VECTOR MULTIPLICATION.
C     THE SUPPLIED MATRIX AMX MUST BE THE TRANSPOSE, TO ENABLE A
C     COLUMN-BY-COLUMN MULTIPLICATION.
C
C********************************************************************
C
      DIMENSION AMX(NR,NC),VEC(NR),RES(NC)
C
      DO 20 I=1,NC
      X=0.E0
      DO 10 J=1,NR
   10 X=X+AMX(J,I)*VEC(J)
   20 RES(I)=X
C
      RETURN
      END
